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ABSTRACT 

We use high resolution simulations of isolated dwarf galaxies to study the physics 
of dark matter cusp-core transformations at the edge of galaxy formation: M200 = 
10^ — 10^ Mq. We work at a resolution (^4 pc minimum cell size; ^250 Mq per particle) 
at which the impact from individual supernovae explosions can be resolved, becoming 
insensitive to even large changes in our numerical ‘sub-grid’ parameters. We find that 
our dwarf galaxies give a remarkable match to the stellar light profile; star formation 
history; metahicity distribution function; and star/gas kinematics of isolated dwarf 
irregular galaxies. Our key result is that dark matter cores of size comparable to the 
stellar half mass radius ri/2 always form if star formation proceeds for long enough. 
Cores fully form in less than 4 Gyrs for the M200 = 10 ^ Mq and ^ 14 Gyrs for the 
10^ Mq dwarf. We provide a convenient two parameter ‘coreNFW’ htting function 
that captures this dark matter core growth as a function of star formation time and 
the projected stellar half mass radius. 

Our results have several implications: (i) we make a strong prediction that if 
AGDM is correct, then ‘pristine’ dark matter cusps will be found either in systems that 
have truncated star formation and/or at radii r > ri/2; (ii) complete core formation 
lowers the projected velocity dispersion at ri/2 by a factor ^ 2 , which is sufficient to 
fully explain the ‘too big to fail problem’; and (hi) cored dwarfs will be much more 
susceptible to tides, leading to a dramatic scouring of the subhalo mass function inside 
galaxies and groups. 

Key words: cosmology: dark matter, galaxies: dwarf, galaxies: haloes, galaxies: 
kinematics and dynamics, methods: numerical 


1 INTRODUCTION 


Dwarf galaxies provide a unique test-bed for galaxy forma¬ 
tion and dark matter physics (e.g. McConnachie|20'T2 ). The 


very smallest - the dwarf spheroidal (dSph) galaxies of the 
Local Group - are resolved in exquisite detail, providing de¬ 


tailed star formation histories (e.g. 

Weisz et al.|2014); orbits 

(e.g. 

Lux et al.||2010); masses (e.g. 

Walker et al.[|2009 Wolf 

et al. 

2010|); and dark matter mass prohles (e.g.jKleyna et al. 

2003 

Goerdt et al.||2006| [Cole et al.||2012| [Battaglia et al. 

2008 

Walker & Penarrubia[[2011 Battaglia et al. 2013). In 


particular, the shallow potential wells of the dSphs are par¬ 
ticularly sensitive to energetic feedback from supernovae. 


providing unique constraints on feedback physics (e.g. Dekel 
Silk||1986| Read et al.||2006b Teyssier et al.|[20T3 ); while 


their high dark matter content makes them natural dark 


matter laboratories (e.g. 

Dalcanton & Hogan 200 1[ Shao 

et al.|2013 

Horiuchi et al. 

2014), and excellent sites to hunt 


E-mail: justin.inglis.read@gmail.com 


for dark matter annihilation or decay (e.g. Lake|19^ [Char- 


2015 


2015 


bonnier et al.||20TT] Malyshev et al.||2014| [Bonnivard et al. 


Drlica-Wagner et al.||20T^ Geringer-Sameth et al. 


It has been known for over a decade that tiny gas rich 


dwarfs favour central dark matter cores over cusps (Flores 


& Primack[ 

[1994[ Moore [1994F a fact that has stood the 

test of tim^ 

(e.g. Kuzio de Naray & Kaufmann 

2011 Oh 

et al.|2011 1 

Hague & Wilkinson|2013). The situation for the 


gas-poor dSphs is more murky, with claims of cores (e.g. 


Battaglia et al. 

2008[ [Walker & Peharrubia[[2011[ [Amorisco 

& Evans! 20111 

Agnello & Evans [2012); cusps (Richardson 

& Fairbairn[[2014 Strigari et al.[[2014[ e.g.); and everything 

in between (e.g. 

Breddels & Helmi|2014 Jardel & Gebhardt 


^ There may be some appreciable scatter, however ( [Simon et al.| 
2005 Adams et al.|2014 Hague fc Wilkinson|2014 ). 
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2013). However, the bri ghtest Milky Way dSph, Fornax, is 
almost certainly corecj^ ( Cole et ^^ 2012 ). 

Such dark matter cores are not predicted by pure 
collisionless cold dark matter structure formation simula¬ 
tions (e.g. Dubinski Carlberg| 199 1| [Navarro et al.|1996b 
Springel et al.|20'08 Stadel et aL|2009| , which has le d many 


to claim that this is evidence for new physic^ 


le.g. 


Moore 


1994 Rocha et al.||2013 Elbert et al.||2015| ). However, dark 


most massive before infall onto the Milky Way. They called 
this the ‘too-big-to-fair problenj^ 

Several authors have noted that if dark matter halos 
are cored rather than cusped then this can naturally allevi¬ 
ate both the too-big-to-fail and missing satellites problems. 
Cored dwarfs have a lower central stellar velocity dispersion 
that - if unmodelled - will cause halo masses to be under¬ 
est imated_([ReaT^r^l][200^ iMadau et al.||2014 Di Cintio 


matter cores can also arise due to energetic feedback from et al.|[^14[ Ogiya & Burkert [2015 ). Such cored dwarfs are 


supernovae, as first suggested by Navarro et al. (1996a). 


Gnedin & Zhao ( 2002 ) showed that the effect for a single 
burst is very small, suggesting that stellar feedback can¬ 
not significantly alter a dark matter cusp. However, |Read| 
& Gilmore (2005) showed that if star formation proceeds in 


repeated bursts then the effect accumulates, gradually grind¬ 
ing a cusp down into a core. This mechanism has now been 
observed in hydrodynamic simulations that resolve the most 
massive sites of star fo rmation (e.g.[Mashchenko et al.|2008| 


Governato et al.|| 2010 | [Teyssier et al.||2013| |Di Cintio et al'l 

2014||Trujillo-GQmez et al.|2015| Ohorbe et al.|2015[ and for 

Pontzen Governato||2014| ), while the physics 


a review see 


of such cusp-core transformations is now well understoo(|^ 


( Pontzen Governato|2012| Pontzen et al.|20T5 ). 

Cusp-core transformations may also be key to solving 
two other small scale puzzles: the missing satellites problem 
(Moore et al. 1999 Klypin et al.|1999' ); and the ‘too-hig-to- 


problem (^ad et al.|2006c Boylan-Kolchin et al.|2011 ). 

The former is a large discrepancy between the number of 
bound dark matter structures predicted in ACDM and the 
number of satellites observed in the Local Group of galaxies. 
The latter is an inconsistency between the abundance and 
mass of satellite galaxies in AGDM and the Local Group. It 
has long been known that the missing satellites problem can¬ 
not be solved by simply pushing stars into the most massive 
dark matter satellites. This results in central stellar velocity 
dispersions that are too high to be consistent with the ob¬ 


served dwarf population (Read et al. 2006c). Boylan-Kolchin 


et al. ( 2011 ) showed that this problem remains even for more 
complex mappings from light to dark - whether placing the 
dwarfs in the most massive halos before reionisation, or the 


^ For Fornax, arguably the most robust constraint comes from its 
globular cluster (GC) system. Its GCs should fall to the centre of 
the dwarf in less than a Hubble time, creating a nucleated dwarf 
that is inconsistent with Fornax’s observed light profile. This can 
be avoided if there is a central constant density core that causes 
dynamical friction to stall ( [Goerdt et al.|2006||Read et al.|2006a| ). 
Such a timing argument has been shown to be remarkably difficult 
to circumvent; it holds even if Fornax is triaxial or affected by 
tides (Gole et al.|2012 Kowalczyk et al.|2013|). 

^ Note that 'warm dark matter' (WDM) models are often invoked 
as an explanation for these cores. However, WDM models cannot 
produce cores of the size required without f ine-tuning (|Strigari| 


i.| 2012 |). 

ith colli- 


et al. 2007 Villaescusa-Navarro &: Dalai 2011 Maccio et al.[ 

^ An alternative mechanism that may act in tandem witti 

sionless heating from supernovae is angular momentum transfer 
from dense infalling lumps (|E1-Zant et al.||2001| |Goerdt et al.| 
[20T0l[C^et al.|2011|[Nipoti Sz Binney|2015| ). Since many cored 
dwarfs are observed to be bulgeless, such dense lumps must then 
be removed from the centre. Stellar feedback can achieve this if 
the dense lumps are purely gaseous (e.g. Nipoti fc Binney|2015| ). 


also much more efficiently tidally stripped, which can signif- 


2006cl [Peharrubia et al.||20101 [Zolotov et al.[[2012[ [Brooks 

et al.[[2013[). 

5 increasing likely that cusp-core 
bursty star formation occur in na- 

While it is becoming 
transformations driven by 

ture (Leaman et al.[[2012 

Teyssier et al.[[2013[ 

Weisz et al. 

2012a Kauffmann 2014 

McQuinn et al. 2015 

), there has 


been quite some debate about the mass scale at which such 


processes become inefficient. Peharrubia et al. (2012) use 


simple energetics arguments to point out that below some 
critical stellar mass, there simply will not be enough inte¬ 
grated supernovae energy to unbind a dark matter cusp. 
Garrison-Kimmel et al. (2013) argue that this mass scale is 


sufficiently large that the cusp-core and too-big-to-fail prob¬ 
lems cannot be solved by stellar feedback. However, |Madau| 
et al. ( 2014| point out that it is easier to unbind cusps at 
high redshift, while Maxwell et al. (2015) show that the as¬ 


sumed radial dark matter profiles matter. They find that 
just a few percent of the supernovae energy is sufficient to 
unbind cusps, even in dwarf galaxies like Fornax (see also 
Amorisco et al. 2014 and Breddels et al. 2015). Such cal¬ 


culations rely on a wide array of assumptions: the mean 
stellar mass to halo mass ratio; the size of the core; the cou¬ 
pling efficiency of the supernovae; the initial mass function 
of stars; the redshift of cusp-core transformations; and the 
dark matter radial profile, amongst others. Differences in 
these assumptions are likely responsible for the wide range 
of results reported in the literature; we will discuss this fur¬ 
ther in 14.11 

Alongside such energy arguments, several groups have 
used hydrodynamical simulations to study cusp-core trans¬ 
formations. These avoid some of the assumptions necessary 
in the above energy arguments, but nonetheless there is 


significant scatter between groups. Di Gintio et al. (2014) 


and Toilet et al. (2016) report a peak in core formation at 
~ 10 ^^ M 0 , with core formation ceasing below this mass 


scale; while Mashchenko et al. (2008), Madau et al. (2014) 


Ohorbe et al. (2015) and Verbeke et al. (2015) find sub- 


slantial cores even in 10^ M© dwarfs. Finally. 


Laporte k, 


Peharrubia (2015) have recently argued that cusps could re¬ 


form as a result of late minor mergers. All of this leaves 
open the question of where, if anywhere, cusp-core transfor¬ 
mations really cease. 

In this paper, we simulate isolated dwarf galaxies over 
the mass range M 200 = 10^ — 10^ M© to model cusp-core 
transformations at the very edge of galaxy formation. Such 
low mass dwarfs have been modelled at high redshift, typi¬ 
cally with a view to understanding the formation of the ‘first 


^ This refers to the fact that it is odd that such massive dwarfs 
have apparently failed to form stars. 
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galaxies’ (e.g. 

Yoshida et al.[[2003 

[Read et al. 

[2006b 

Wise 

& Abel||2007| 

Holey et al.[[2009[ [J 

ohnson et a 

112009 

Wise 


neously resolved stellar feedback and followed the evolution 
of such dwarfs over a Hubble time. This is the goal of this 
present work. All simulations presented in this paper reach 
a minimum cell size of ~ 4 pc, allowing us to resolve the 
cooling occurring in the majority of individual supernovae 
remnants. This allows us to correctly capture the momen¬ 
tum injection into the interstellar medium (ISM) generated 
by expanding bubbles of shock heated gas. As we will show, 
this makes us insensitive to the details of our ‘sub-grid’ nu¬ 
merical parameters (for a detailed discussion on this, see e.g. 
[Kirnm et al. 12015) for earlier work on the utility of resolving 


the interstellar medium (ISM) see e.g. Kravtsov 2003 [Saitoh 
et al.|[2008 Hopkins et al.|[^011| 2013 Agertz et al.||2013 ). 
For this reason, we consider the simulations presented here 
as predietive; they allow us to calculate how dark matter 
density profiles evolve in the face of stellar feedback, from 
first principles. Although we miss some potentially impor¬ 
tant physics - mergers; tides; ram pressure stripping and 
photoionisation - we show that our simulated dwarfs give a 
remarkable match to all known observational data for real 
isolated dwarfs in the Universe. This suggests that our sim¬ 
ulations already capture the essential physics - at least for 
isolated dwarfs. 

This paper is organised as follows. In ^ we describe 
the simulation suite. In E we present our key results. We 
give an overview of the simulations (^3.1); compare them 
with a wide array of data for two isolated dwarfs - Leo T 


and Aquarius (J3.2); and present our results for dark matter 
cusp-core transformations (J3.3). We show that at the reso¬ 
lution we adopt here, our results are insensitive to our choice 


of ‘sub-grid’ numerical parameters (f 3.3.1) and only weakly 
sensitive to our choice of initial conditions (f 3.3.2). We de¬ 
rive a convenient COReNFW fitting formula that captures 
the evolution of the dark matter cusp as a function of the 
projected stellar half mass radius, and the star formation 
time (f3.3.3). Finally, in ^ and we discuss our results 
and their implications for cosmology, and present our key 
conclusions. 


2 THE SIMULATIONS 
2.1 Initial conditions 


We set up the dark matter halos following [Read et al.j 
(2006c). The particles were populated using accept/reject 
from an analytic density profile; their velocities were drawn 
from a numerically calculated distribution function, assum¬ 
ing an isotropic velocity dispersion tensor. We assume a 
Navarro et al. ( 1996b| (hereafter NFW) form: 


PNPw(r) = po (^) (l + ^) (1) 

where the central density po and scale length are given 
by: 

Po — Pcrit^C Qcl^ 5 — '^20o/c (2) 


and 


9c = 


1 

log(l + c) 


c 

1+c 


^200 = 


3,^ 1 

yM2oo—r- 

4 TrApcrit 


1/3 


(3) 

(4) 


where c is the dimensionless eoneentration parameter] A = 
200 is the over-density parameter; pcHt = 136.05 M© kpc“^ 
is the critical density of the Universe at redshift z = 0; r 2 oo 
is the ‘virial’ radius at which the mean enclosed density is 
A X pcrit; and M200 is the ‘virial’ mass - the mass within 
^200. 

The enclosed mass for the NFW profile is given by: 


Mnfw(^) = M2oogc 


In 




(5) 


which logarithmically diverges as r ^ 00. We truncate the 
mass distribution at r = 3 r2oo • 

We consider four halos with virial masses M200 = 
10^, 10®, 5 X 10® and 10^ M©; these have dark matter virial 
masses of M2 oo,dm = M 2oo/(l — fb), where fb = 0.15 is the 


Universal baryon fraction (e.g. [Planck Collaboration et "al 


2014). Their concentration parameters were selected to be 


the cosmic mean value taken from cosmological simulations, 
following Maccio et al.[ ( [2007[ ) . 

These halos were then filled with the Universal baryon 
fraction in gas {fb = 0.15), set up either as an NFW pro¬ 
file in hydrostatic equilibrium or as a constant density slab 
in hydrostatic equilibrium (in order to explore our sensi¬ 
tivity to this choice). The gas was given a seed metallicity 
of Zgas = 1 O~®Z0 , representing Pop HI enrichment (e.g. 


Nakamura Umemura[[^001[ [Holey et al.[[20091 [Karlsson] 


et al. 


2013). Experiments with the initial Zgas in the range 
< ^gas/^o < 10~^ show no significant differences af- 
2 Gyr, as self-enrichment after the first burst 


10 “ 

ter the first 

of star formation becomes the dominant source of metals. 
Note that such differences may matter more in the presence 
of the UV background; we will consider this in a forthcoming 
paper. 

We added angular momentum to the gas assuming 
a specific angular momentum profile as in [Bullock et al.[ 

(I2OOTI): 


( 6 ) 

M 2 OO 

where the peak specific angular momentum jmax is set such 
that the total halo angular momentum is given by: 


/ O// 


(7) 

( 8 ) 


J 200 = 47r / j{r)pT^FW(r) 

Jo 

= \'yj2GMiooR2oo 

where A' is the spin parameter. We assume the cosmic mean 
value A' = 0.035, except for one simulation where we explore 
the effect of higher A' = 0.07. 

The default particle number was chosen such that each 
dark matter particle had mass ttidm = 250 M©, similar to 


the stellar masses (see ^2.2); one simulation was run at ten 
times this resolution to test for numerical convergence. We 
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explicitly checked that our initial conditions remain in equi¬ 
librium for a Hubble time (see Appendix [A|. All simulation 
parameters and our simulation labelling system are sum¬ 
marised in Table [T] 

The above choice of ‘cooling halo’ initial conditions has 
been widely used in the literature to study galaxy forma¬ 
tion (e.g. Kaufmann et al.|20Q6 Teyssier et al.|20T3 Hobbs 


2008). Each star particle then looses mass due to supernovae 


et al.|20l5 ). It has the advantage that it is well-defined and 

therefore useful for controlled numerical experiments. It is 
also likely to be realistic for tiny dwarf galaxies that un¬ 


dergo few mergers over most of their lifetimes (e.g. Verbeke 
et al.|[20T5 ). However, such initial conditions have the dis¬ 
advantage that they start out with a dense central concen¬ 
tration of gas that leads unavoidably to an initial starburst 


(see e.g. discussion in Hobbs et al. 2015). To explore the 
effect of this on our results, we ran an additional simula¬ 
tion: M5e8c25_2e6_rhocon. This is identical to our fiducial 
M200 = 5 X 10® M© simulation (M5e8c25_2e6), but start¬ 
ing out instead with a constant gas density out to the virial 
radius. We discuss the results of this simulation in §3.3.2[ 


2.2 Sub-grid physics 

Our suite of simulations were carried out using the Adaptive 


Mesh Refinement (AMR) code RAMSES (Teyssier 2002). 
The fluid dynamics is calculated using a second-order unsplit 
Godunov method, while the collisionless dynamics of stellar 
and dark matter particles is evolved using the particle-mesh 
technique ( Hockney fc Eastwood|[l981 ), with gravitational 
accelerations computed from the gravitational potential on 
the mesh. The potential is calculated by solving the Poisson 


equation using the multi-grid method (Guillet & Teyssier 


2011) for all refinement levels. The equation of state of the 


fluid is that of an ideal mono-atomic gas with an adiabatic 
index 7 = 5/3. 

Our refinement strategy is based on a quasi-Lagrangian 
approach; each cell is refined if it contains more than 8 dark 
matter particles, or if its baryonic mass (including gas and 
star particle mass) exceeds 8 x rrires, where rrires = 60 M© is 
the adopted mass resolution of the simulation. This allows 
the local force softening to closely match the local mean 
interparticle separation, which suppress discreteness effects 
(e.g. Romeo et al. 2008). Refinement is performed recur¬ 


sively, on a cell-by-cell basis, until the adopted maximum 
allowed level of refinement is reached. In our current setup, 
the finest grid cell size is Ax 4 pc. 


2.2.1 Star formation 

The adopted star formation and feedback physics are pre¬ 


sented in detail in Agertz et al. (2013) and Agertz & 


Kravtsov (2015a). Briefly, we adopt a local star formation 


rate using a Schmidt relation 


= eff —, p> 


(9) 


where pg is the gas density in a cell, tff = f32Gpg the 
local free-fall time of the gas, and Cff is the star formation 
efficiency per free-fall time. Eq.|^ is sampled via a Poisson 
process with star particles masses being multiples of the 


explosions and stellar winds, see |2.2.2| 

Star formation is known to correlate well with the pres¬ 
ence of molecular gas (e.g. Bigiel et al.|2008 ), i.e. H2 traced 


observationally by GO. We note that star formation from H2 
in low metallicity environments (Zgas ^ 10“^ Zq) is poorly 
understood, and difficult to model, due to uncertainties in 
dissociating UV radiation; self-shielding and shielding of H2 
by dust ( Krumholz et al.|[2QQ9 Gnedin et al.||2Q09 Gnedin 


fc Kravtsov|201 1 ); uncertain dust-to-gas ratios (Eisher et al. 

2014); and the physics of line overlap (Gnedin & Draine 


2014). Eor these reasons, we adopt here a high star forma¬ 
tion density threshold, = 300mHcm“®, ensuring that 
star formation proceeds only in gas dense enough to be 
mainly molecular gas. We explore the effect of raising this 
to p* = 10® mu cm“® in one run. 

In our fiducial models, we adopt a star formation effi¬ 
ciency per free fall time eff = 0.1 motivated by observations 
of local GMGs ( Lada et al.|2010 Murray [2011 Evans et al 


2014), as well as recent numerical work ( Hopkins et al.|2014 
Agertz Kravtsov|2015a ) demonstrating how the interplay 


between local efficient star formation and stellar feedback 
reproduces the low global star formation efficiency inferred 
from the Kennicutt-Schmidt relation (e.g. Bigiel et al.|2008 ) 
and global galaxy scaling relations. In one run, we explore 
the effect of lowering this by a factor of ten. 

2.2.2 Stellar Feedback 


We adopt the stellar feedback model described in |Agertz| 
et al. (2013). Briefly, each formed stellar particle is treated 


as a single-age stellar population with a Ghabrier (2003) ini¬ 


tial mass function (IME). We account for injection of energy, 
momentum, mass and heavy elements over time via super¬ 
nova type H (SNH) and type la (SNIa) explosions, stellar 
winds and radiation pressure (allowing for both single scat¬ 
tering and multiple scattering events on dust) on the sur¬ 
rounding gas. Each mechanism depends on the stellar age, 
mass and gas/stellar metallicity, calibrated on the stellar 
evolution code STARBURST99 ( Leitherer et al.|1999 ). Eeed- 
back occurs continuously at the appropriate times when the 
various feedback process are known to operate, taking into 
account the lifetime of stars of different masses in a stellar 
population. To track the lifetimes of stars within the popu¬ 
lation we adopt the metallicity dependent age-mass relation 
of Raiteri et al. (1996). 


Momentum from stellar winds, radiation pressure, and 
SNe blastwaves is added to the 26 nearest cells surrounding 
a parent cell of the stellar particle. Thermal energy from 
shocked SNe and stellar wind ejecta is injected directly into 
the parent cell. In the case of SNe we model these as dis¬ 
crete events (i.e. we sample the stellar IME discretely) and 
add 10^^ ergs per SN to the gas thermal energy. As the 
average gas and stellar metallicities in our suite of simu¬ 


lated dwarf galaxies are low (Zgas < 0.1, see §3.2.3), we 
expect the impact from stellar winds and radiation pressure 
on dust to be relatively small compared to that of the SNe. 
Heavy elements (metals) injected by supernovae and winds 
are advected as a passive scalar and are incorporated self- 
consistent ly in the cooling and heating routine. The code ac¬ 
counts for metallicity dependent cooling by using tabulated 


sampling mass m* = 300 M© (see e.g. Dubois & Teyssier cooling functions of Sutherland & Dopita (1993) for gas tern- 
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Label 

M200,DM 

[Mq] 

C 

A^(< ''’200) 

A' 

fb 

[2.33 X 10® Mq kpc km/s] 

J 200 

[kpc km/s] 

imax 

Description 

M7c36_4e4 

10 ^ (1 - /(,) 

36.73 

4 X lO'^ 

0.035 

0.15 

29.47 

1.35 

Fiducial 

M8c28_4e5 

10* (1 - A) 

28.57 

4 X 10® 

0.035 

0.15 

1367.88 

6.28 

Fiducial 

M8c28_4e5_e001 

10* (1 - A) 

28.57 

4 X 10® 

0.035 

0.15 

1367.88 

6.28 

Cff = 1% 

M8c28_4e5_pST 

io» (1 - A) 

28.57 

4 X 10® 

0.035 

0.15 

1367.88 

6.28 

Sedov-Taylor momentum (pst) 

M8c28_4e6_rle3 

io» (1 - A) 

28.57 

4 X 10® 

0.035 

0.15 

1367.88 

6.28 

Higher DM resolution; — lO^mncm"^ 

M5e8c25_2e6 

5 X IQS (1 - A) 

24.93 

2 X 10® 

0.035 

0.15 

19996.3 

18.38 

Fiducial 

M5e8c25_2e6_rhocon 

5 X 108 (1 - A) 

24.93 

2 X 10® 

0.035 

0.15 

19996.3 

18.38 

Constant initial gas density over 0 < r < r2oo 

M9c22_4e6 

10 '' (1 - A) 

22.23 

4 X 10® 

0.035 

0.15 

63481.41 

29.17 

Fiducial 

M9c22_4e6_lam007 

10 '' (1 - A) 

22.23 

4 X 10® 

0.07 

0.15 

126962.81 

58.34 

Double fiducial halo spin (2 x A') 


Table 1. The simulation suite. The columns show from left to right: the simulation label; the initial NFW halo parameters (see equation 
[T] >; the dark matter particle number N within the virial radius r2oo; the spin parameter A'; the baryon fraction the total angular 
momentum J200; the maximum specific angular momentum jmax; and a description of the simulation. For further details, see 0 


peratures 10^ — 10^'^ K, and metal f i ne-str ucture cooling be¬ 
low 10^ K as in Rosen & Bregman (1995). 

At the current numerical resolution {Ax ~ 4pc), and 
for the density field and typical metallicities in these sim¬ 
ulations, we are likely to resolve the impact of most SNe 
explosions. For a point like explosion, the cooling radius, 
i.e. the bubble radius for which radiative losses are ex¬ 
pected to be important for each discrete SN event, scales 
as rcooi ~ SOuq °'^^(Zgas/^o pc for a supernova 

explosion with energy F^sn = 10 ^^ erg and ambient gas den- 
sity no (in units of cm“^; e.g. Cioffi et al.||l988 Thornton 


|et al.|p!9^ Kim Ostriker||2015 ). Kim & Ostriker (2015) 


(and see also Martizzi et al.|20T5 Gatto et al.|2015 Simpson 

|et al.|20l5 ) demonstrated that if rcooi is resolved by at least 
three grid cells (rcooi ^ 3Ax), the correct amount of mo¬ 
mentum generated in the Sedov-Taylor phase is recovered 
before the bubble energy is radiated away. 

For typical ISM gas densities, thi s momentum, ng-r ^ 
2.6 X 10^ F^ 5 i^^^nQ M© kms“^ (e.g. Blondin et al.|l998 ), 
can be > 10 times greater than the initial ejecta momen¬ 
tum, which is why it is important to capture this process. 
For Ax = 4 pc, and Zgas = 10~^ Zq, we resolve single ex¬ 
plosions up to ambient densities of ~ 50cm“^. Type II SNe 
occur for almost ~ 40 Myr (roughly the age of an 8 M© star) 
in our models, and during this time pre-SNe feedback and 
gas turbulence preprocesses the gas, allowing for explosions 
to be resolved. In §3.3. 1| we demonstrate that this is the 
case, and that our conclusions are insensitive to wether we 
explicitly inject psT, rather than allowing it to be generated 
“above the grid” via individual expanding SNe bubbles. 


3 RESULTS 
3.1 Overview 

In Figure ^ we give an overview of our simulation results 
over the mass range 10® — 10^ M©. From top to bottom, 
the rows show the projected gas density pgas’, temperature 
^gas; and the stellar density p* for an edge on view of the 
galaxies. The columns show the 10® M©, 5 x 10® M© and 
10^ M© simulation, respectively. Notice that very few stars 
form in the 10® simulation - just 86 star particles - which 
is too few to make a plot of p*. By contrast, the 10® M© 
simulation forms over 15,000 star particles, and overlapping 
SNe explosions drive large hot bubbles in the disc as can be 
seen in panels (c) and (f). 

We do not show results for the 10^ M© simulation since 


this forms just 1-2 star particles after which star forma¬ 
tion ceases. This agrees with earlier works that suggest 
that stellar feedback sets an edge to ga laxy formation at 
M200 ~ 10^ M© (e.g. Read et al. 2006b). However, unlike 
these previous works, here galactic winds are not modelled 
in a sub-grid fashion, but emerge self-consistently by resolv¬ 
ing individual interacting SNe explosions. 

We note that the precise edge of galaxy formation will 
be very sensitive to the details of reionisation; epoch of gas 
infall; star formation; cooling physics; and ‘population IIF 
star formation/evolution that we do not properly resolve or 


model here (e.g. Efstathiou 1992[ Ricotti et al. 2008 Johnson 
et al.| 200^ [Boley et al.||2009 [Wise et al. 2012). For this 


reason, while we are confident that galaxies at 10^ M© will 
contain very few stars, we defer further exploration of this 
mass scale to future work where we will consider such physics 
more carefully. 


3.2 Observational properties 


Before discussing the evolution of the dark matter distribu¬ 
tion in our simulations, we first compare them with a host of 
observational data for two isolated field dwarfs: Leo T and 
Aquarius (also called DDO 210). These are two of the small¬ 
est galaxies in the Universe with on-going star formation 


dlrwin et al.||2007 Ryan-Weber et ^|2Q08| van den Bergh 


1959|). At 4 17^19 kpc (|Simon Geha||2Q07| and ~ 1 Mpc 
( |Lee et ^[1999 ) from any large spiral, they are extremely 
isolated. Both have measured star formation histories (Weisz 


et ah 2012b Gole et ah 2014); gas kinematics (Begum & 

Ghengalur||2004 Ryan-Weber et ah 2008 Oh et ah 

2015); 

stellar kinematics (Simon & Geha||2007 Kirby et ah 

2014); 

photometry (Martin et aL||2008| McGonnachie et ah 

2006); 


and stellar metallicity distribution functions (Weisz et al. 
2012b Kirby et al.||2013| Gole et ^|2014 ). For this reason, 
they are the natural first place to look for data comparisons 
with the simulations we present here. The results from the 
M200 = 5 X 10® M© and 10® M© simulations (see Table[^ 
are compared to observational data for Leo T and Aquar¬ 
ius in Figures and The panels show the stellar surface 
density (a); star formation rate (b); stellar metallicity dis¬ 
tribution function (c); gas rotation curve (d); gas velocity 
dispersion (e); and the projected stellar velocity dispersion 
(f). We discuss these in turn, next. 
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- 2-10 1 2 - 2-10 1 2 
x/kpc x/kpc 


Figure 1. Overview of the simulation results for our fiducial simulations: (a-c) projected gas density Pgas; (d-f) gas temperature Tgas] 
and (g-h) stellar density p*, for an edge on view of the galaxies. The columns show the M 200 = 10® Mq, 5 x 10® Mq and 10^ Mq fiducial 
simulations (M8c28_4e5; M5e8c25_2e6; M9c22_4e6; see Table [l| , respectively. Notice that very few stars form in the 10® simulation - 
just 86 star particles - which is too few to make a plot of p*. By contrast, the 10® Mq simulation forms over 15,000 star particles, and 
overlapping SNe explosions drive large hot bubbles in the disc as can be seen in panels (c) and (f). We do not show results for the 10^ Mq 
simulation since this forms just 1-2 star particles after which star formation ceases. 


3.2.1 Projected light profiles 


The vertical green dashed lines in Figures]^ a) and|^a) mark 
the projected 2D half mass radius of the simulations; the red 
bands mark the similar confidence intervals for LeoT (|Mar- 
tin et al. 2008) and Aquarius (McConnachie et al. 2006), 


respectively. The red and magenta data points show data 
from resolved star counts for LeoT ( Martin et al.|20^ ) and 
Aquarius (Martin et al. 2008); we renormalise these to match 


our simulated surface mass density at the half stellar mass 
radius. The brown data points in Figurej^a) show the stellar 
surface mass density for Aquarius derived from integrated 
light ( Zhang et al.|2012 ). This is also renormalised to match 


our simulated profile at the projected half stellar mass ra- 

diu0 

The data in McConnachie et al. (2006) for Aquarius are 


® If we do not renormalise, then we find that the observed stellar 
surface density profile is lower than our fiducial M 200 = 10® Mq 
simulation by a factor ~ 5 — 7. However, there is at least a factor 
~ 2 uncertainty in the stellar mass of Aquarius ( [Kirby et al.|201^ 
[Zhang et al.|20i^ , while our simulations overproduce stars in the 
first 2 —4 Gyrs due to our idealised ‘cooling halo’ initial conditions 
(see discussion in p.l[ and ^3.2.2[ ). If we exclude stars from the 
transient start-up phase, and take the larger [Kirby et ^ ( [2013] ) 
stellar mass for Aquarius, then our 10® Mq fiducial simulated 
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sufficiently good that the surface brightness could be split 
by age into old and young stars (< GOMyrs; magenta data 
points in Figure]^ a)). Performing a similar cut on our simu¬ 
lations, we qualitatively recover the observed result that the 
young stars (red lines) are more concentrated than the full 
distribution (blue), with a steeper surface brightness fall-off. 
This occurs because the stars are collisionlessly heated by 


stellar feedback, similarly to the dark matter (see J3.3 and 


.g.|Read Gilmore|2005||Leaman et al.|2012||Teyssier et al.| 


2013||E1-Badry et al.|201^ Gonzalez-Samaniego et al.|2015| ). 

Over time, this leads to the older population becoming hot¬ 
ter and more extended than the younger stellar population. 
Such heating occurs very rapidly, as shown by the red lines 
on panel Figure]^ a). These mark the stellar surface mass 
density for two age cuts at < 100 and < 200Myrs. Notice 
that the stars < 200 Myrs old are already substantially more 
extended than those < 100 Myrs old. 

While our M 200 = 10^ M© fiducial simulation gives a 
remarkable match to the real data for Aquarius, our M 200 = 
5x10® M© fiducial simulation does not rise as steeply as the 
data for Leo T (compare the red data points and solid blue 
line in Figure [^a)). Interestingly, however, if we incline our 
simulation at 20° (i.e. near face-on), then it does give an 
excellent match (blue dashed line). We will return to this 
when discussing LeoT’s rotation curve and projected stellar 
velocity dispersion, below. 


3.2.2 Star formation histories 

The star formation histories are given in Figures and 
panels (b). These show the star formation rate (SFR) as a 
function of time t, where t = 0 is the beginning of the Uni¬ 
verse. The blue lines show the raw simulation results; the red 
and purple show data for Leo T and Aquarius, respectively 
(taken from Weisz et al. 2012b and Cole et al. 2014); the 
green data points show the simulation results binned simi¬ 
larly to the observational data for Leo T. The bursty nature 
of the star formation in the simulations is evident; this is 
what is ultimately responsible for the dark matter cusp-core 


transformations that we will discuss in f3.3 (c.f. Read & 
Gilmore[|2005| Pontzen Governato||2012 ). However, when 


binned at the age resolution of observed star formation his¬ 
tories (that come from fits to the colour magnitude diagram 
of stars), the bursts are much harder to spot. At this age res¬ 
olution, the simulations - the like observational data - look 
very smooth, with a near-continuous star formation history 
over cosmic time. 

Our 5 X 10® M© simulation gives a remarkable match 
to the data for Leo T, while the 10^ M© simulation gives a 
much better match to Aquarius. This suggests that the star 
formation history is very sensitive to the halo mass. Indeed, 
we may imagine using the star formation history, calibrated 
on these models, as a sensitive indicator of M 200 in low mass 
systems. This is potentially very powerful as it can work as 
a mass indicator even for dwarfs that have had their star 
formation truncated on in-fall to a galaxy or group. We will 
discuss such ideas further in forthcoming publications. 

While the mean star formation rate agrees very well 


dwarf gives an excellent match to the stellar surface mass density 
profile for Aquarius in both its normalisation and its shape. 


between the simulations and the data, there are differences 
in the fine structure. Looking at Figures]^ andpanel (b), 
we can see that both simulations overproduce stars over the 
first ~ 2 — 4 Gyrs with respect to Leo T and Aquarius (com¬ 
pare the green data points with the red/purple data points). 
This owes to our idealised initial conditions. By starting out 
with a fully-formed dark matter halo filled with dense gas, 
an initial starburst is inevitable. We discuss the sensitivity 
of our results to our choice of initial gas density distribution 
in j jO:2l 

In real galaxies, early star formation will be suppressed 
both by reionisation and mergers, neither of which are mod¬ 
elled here. The former pre-heats gas, suppressing cooling 
and star formation; the latter builds galaxies out of smaller 
sub-units that form stars less efficiently. Indeed, Aquarius 
shows statistically significant evidence for an enhancement 
in star formation ~ 6 Gyrs ago that could owe to some late 
gas accretion (possibly due to a merger; see the discussion 
Gole et al. 2014). We discuss these issues further in 0 


but note here that such differences between the theoretical 
models and the data are much smaller than the difference 
in mean star formation rate between Aquarius and Leo T 
(or equivalently, between our M 200 = 5 x 10® M© fiducial 
simulation and our M 200 = 10^ M© simulation). 

3.2.3 Stellar metallicity distribution funetions 

In Figures [^c) and[^c), we compare our simulated stellar 
metallicity distribution functions (MDFs; blue histograms) 
with data for Leo T and Aquarius (red bands). For the data, 
we use the mean metal content [M/H] as determined from 
fits to the colour magnitude diagram. This measure is clos¬ 
est to our simulated MDFs that are all scaled relative to 
solar abundance. The spectroscopic measures of [Fe/H] are, 
however, in both cases very similar and do not affect the 
conclusions we draw here (Kirby et al. 2013| ). The width 
of the red bands shows the measured variance of [M/H]; in 
both cases, this is also in excellent agreement with the sim¬ 
ulations. 

The above results indicate that galactic winds have reg¬ 
ulated the gas metal reservoir in a realistic manner. To bet¬ 
ter understand the role of metal rich outflows, we calculate 
effective yield, defined as: 


2/eff — 


ln(l//g: 


( 10 ) 


where /gas = Mgas/(Mgas + M*) is the fraction of baryons 
in the gas phase. The effective yield has been widely used 
as a diagnostic of the evolution of the baryonic component 
of galaxies, and more specifically as a test of the validity 


of the closed-box approximation (Pagel & Patchett 1975 


Edmunds 1990). Observationally, the effective yield is known 
to decrease with galactic mass ( Tremonti et al.||20Q4| , with 
a sharp de cline around th e mass of dwarf galaxies (urot ^ 
100kms“^, Garnett|2002 ). 

Under a “closed-box” assumption (no inflow or outflow 
of material), the effective yield is always equal to the true 
yield ytme- This is typically defined for a single stellar pop¬ 
ulation as the mass in newly synthesised metals returned to 
the ISM, normalised to the stellar mass of this population 
in stellar remnants and long-lived stars. Eor our feedback 
prescription, and choice of IME, ytme = 0.022. 
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Following Tassis et al. (2008), (see also Brooks et al. 

2007 1 Agertz & Kravtsov||2015a[ 

), we calculate the observed 
, where we consider only the 

effective yield using Equation] 10| 


cold gas (T ^10^ K) and the metal content within the stellar 
extent, defined as the radius that includes 90% of the total 
stellar mass. The effective yields for the LeoT and Aquarius 
simulations are ^eff ~ 4.8 x 10“^ and 5.5 x 10“^, respec¬ 
tively. These values are lower than the true yield, indicating 
that the gas and stellar metallicities (Figures and are 
controlled by supernovae driven outflows. 


3.2.4 Gas kinematics 

Figures l^d) and|^c) show the rotation velocity of the gas 
V(f),gas (magenta lines). We plot results only for gas with tem¬ 
perature T < 10^ K to mimic HI gas observations. For the 
M200 = 5 X 10® Mq simulation (Figure]^, we plot the rota¬ 
tion curve also at two different inclination angles, 20° and 
60° (dashed lines, as marked; 0° corresponds to a face on 
galaxy). For the 10® Mq simulation (Figure]^, we show the 
time dependance of r’</),gas- The translucent magenta lines 
show r’(^,gas at times [0,0.5,1,1.5] Gyrs ago. As can be seen, 
there is appreciable scatter in r’(/>,gas as a function of time. We 
hnd that this scatter is slightly larger for the high spin ver¬ 
sion of this simulation, M9c22_4e6_lam007. However, since 
the results for M9c22_4e6_lam007 are otherwise very simi¬ 
lar to M9c22_4e6, we omit it for brevity. Also marked on 
the plots are the circular speed curves for the stars (blue); 
dark matter (green); gas (red); and all mass (cyan). For the 
10® M© simulation, we overplot the asymmetric drift cor¬ 
rected rotation curve data for Aquarius (red data points), 
taken from Oh et ah] ( |2015| ) . Panel (e) shows that the gas 
velocity dispersions (blue), as compared to data (red). 

For Leo T (Figure]^, no rotation curve has been re¬ 
ported, but the gas velocity dispersion of the 5x10® M© 
simulation is in excellent agreement with reported measure¬ 
ments ( Ryan-Weber et ar]|2008 red dots). Interestingly, if 
Leo T is inclined at ~ 20° - as would give the best ht to the 
photometric light prohle (see [ 3.2. 1| - then the gas rotation 
would indeed be very difficult to detect, with a peak non¬ 
inclination corrected r’(^,gas of just under ~ 5km/s (magenta 
dashed line). When aligned edge-on, the gas rotational ve¬ 
locity (magenta solid line) gives a reasonable match to the 
circular speed curve (cyan line), with almost all of the mass 
at all radii coming from the dark matter (solid green line). 
The gas and star contributions (red and blue lines) are com¬ 
parable, similarly to what has been reported recently for low 
mass isolated dwarf irregulars in the LittleTHINGS survey 


(Oh et al. 2015). We will compare our models to rotation 


curves in more detail in a companion paper. 

Our 10® M© simulation gives a good qualitative match 
to the observed rotation curve and gas dispersion of Aquar¬ 
ius (see panels (d) and (e) in Figure |^, if aligned edge on. 
(For this reason, we do not consider alternative inclination 
angles for our 10® M© simulation.) However, the amplitude 
of rotation is somewhat higher for Aquarius, suggesting that 
it is either slightly more massive than M 200 = 10® M©, 
and/or has higher concentration. 


3.2.5 Projected stellar velocity dispersions 


Finally, we consider the observed projected stellar velocity 
dispersions. For both Leo T and Aquarius, the latest data 
amount to a single measurement at ~ the projected stellar 
half light radius. This is shown in Figures and panels 
(f). The blue lines show the simulation results aligned edge 
on; the dashed aligned at 20° and 60°, as marked. Notice 
that for Leo T, the 20° inclination gives a reasonable match 
to the data, but is ~ 2km/s lower. The red line shows a ht 
to the simulation data using the Jeans equation, where we 
account also for the fact that the stars rotate signihcantly 
in the x — y plane. Specihcally, we solve the radial Jeans 
equation (e.g. Binney Tremaine||2008 ): 


C^LOS — 


S.(fl) 


/:( 




v^{r)(jl{r)r 

- R 2 


dr 


( 11 ) 


where E*(R) is the surface brightness density of the stars; 
(r) is the spherically averaged stellar density; 


13 = 1- 


cTe 5-crl 

2a^ 


( 12 ) 


is the velocity anisotropy; and (Jr, 0,0 are the radial, 0 and cj 
velocity dispersions, respectively. The radial dispersion re¬ 
lates to the enclosed mass M(r) as: 

^2 ^ /W r 

where: 


/(r) = exp (^-2 J (14) 


and: 


= 1 - 


2 I 2 I —2 

cre5-cr^P v^ 
2 cr^ 


(15) 


where (3' accounts for the mean streaming motion in the 
plane This is an important contributor since our sim¬ 
ulated stellar distributions are signihcantly rotating. This 
rotation explains why the face-on projected velocity disper¬ 
sions are rather low (Figures [^f) and|^f) solid blue lines). 

We use our perfect knowledge of z/*, S*, /3', /3 and M(r) 
to determine a ‘Jeans equation’ ctlos, by numerically solv¬ 
ing equation m this is marked by the red lines on Figures 
l^f) and[^f). For the 5 x 10® M© simulation, this gives an 
excellent match to the face-on 20° simulation data (blue 
dashed line). For the 10® M© simulation, it is slightly off at 
large radii, most likely owing to the fact that the system is 
not quite in a steady-state (a key assumption of the Jeans 
analysis). The agreement is nonetheless good. 

Armed with our Jeans analysis, we can now determine 
what the stellar velocity dispersion would look like in the 
absence of rotation or if we ‘undo’ the cusp-core transforma¬ 
tion. The green lines show the effect of switching off rotation; 
the cyan of enforcing the initial NFW prohle, all other things 
being equal. Removing rotation boosts the projected veloc¬ 
ity dispersion, particularly at large radii. This may be closer 
to the situation in real dwarfs if they undergo signihcant late 
mergers. We do not model such mergers here, but they will 
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Figure 2. A comparison of our M 200 = 5 x 10® Mq fiducial simulation (M5e8c25_2e6) with a host of observational data for the Leo T 
isolated dwarf irregular galaxy. The panels show: (a) The stellar surface density E*. The simulation results for an edge-on projection are 
marked by the blue lines; the dashed lines show the same inclined at 60° and 20°, as marked (where 0° corresponds to a face-on view). 
The vertical green dashed line marks the projected stellar half mass radius (and similarly for all other panels). The red data points show 
the photometric data for Leo T, taken from [Martin et ^ ( |2008] >; the red band marks the observed projected half mass radius, including 
uncertainties, (b) The star formation rate (SFR) as a function of time t (where t = 0 is the beginning of the Universe). The blue line 
shows the full simulation results; the green re-binned to match the Leo T data. Overplotted are measured star formation histories for Leo 
T and Aquarius, taken from jWeisz et al.] ( |2Q12b| >; [Cole et al.j ( |2014| ). (c) Stellar metallicity distribution function, using the scaled Solar 
mean metallicity [M/H]. The histogram shows the simulation results; the mean is marked by the vertical dashed line. The data for Leo 
T are marked by the red band; the width marks the measured dispersion (data taken from jWeisz et al.|2012b| >. (d) Gas rotation velocity 
'^(/),gas (magenta line). The magenta dashed lines mark two different inclination angles (20°; 60°, as marked). The blue, red, green and 
cyan lines mark the circular speed curves (measured directly from the gravitational potential) for the stars, gas, dark matter, and all 
matter, respectively, (e) Gas velocity dispersion (blue). Two data points for Leo T (reported without errors) are overplotted, taken from 
[Ryan- Weber et ^ ( |2008| . (f) Projected stellar velocity dispersion. The simulation data are shown in blue; the dashed lines give results 
for two different inclination angles (20°; 60°, as marked). The red line shows a Jeans equation estimate of ctlos that accounts for the 
significant stellar rotation (see text for details). The green line shows ctlos with this rotation removed. The cyan line shows ctlos if we 
restore the initial dark matter cusp, all other things being equal. Notice that the green line produces a dispersion profile similar to that 
seen if Leo T is nearly edge-on, while reinstating the dark matter cusp (cyan line) boosts the central ctlqs by nearly a factor of two 
(compare the blue and cyan lines). The red data point shows the measured ctlos for Leo T ( [Simon Sz Geha|2007| >. 


act to randomise the rotation, lowering and moving us 
more towards the green lines. Switching to a cuspy NFW 
profile mostly alters the central velocity dispersions inside 
~ '^ 1 / 2 - The effect, however, is dramatic - nearly doubling 
ctlos in both cases (compare the blue and cyan solid lines in 
Figures!^ andpanels (f)). This highlights that using stel¬ 
lar velocity dispersions alone to estimate halo masses could 
lead to highly biased conclusions, if assuming an NFW dark 
matter mass prohle. We discuss this further in Q 

We stress that neither of our simulations has been fit to 
the data. Nor has any ‘hne tuning’ of our simulation param¬ 


eters been performed. Our only freedoms are in the inclina¬ 
tion angle of our simulations; and which galaxy we choose 
to compare our simulations to. For this reason, it is strik¬ 
ing that the agreement of both simulations with Leo T and 
Aquarius, respectively, is so good. It suggests that we have 
captured the essential physics relevant for modelling isolated 
dwarfs in the held. We discuss this further in Q 
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Figure 3. A comparison of the M 200 = 10^ Mq fiducial dwarf (M9c22_4e6) with a host of observational data for the Aquarius dwarf 
irregular galaxy. The panels are as in Figure]^ with two exceptions. In panel (a), we show results also for the young stars (< 100 and 
< 200Myrs old, as marked in red), and similarly for the Aquarius data (< OOMyrs old; magenta data points). In panel (d), we show 
results for at different times [0, 0.5,1,1.5] Gyrs ago (magenta lines). The scatter is caused by the bursty star formation that blows 

large HI bubbles through the disc. Data for the asymmetric drift corrected Aqarius rotation curve (red data points) are overplotted 
(taken from [Oh et al.|2015|>. The data for Aquarius are compiled from: [McConnachie et al.| ( |2006| >; [Zhang et al.| ( |2012| > (ZH12); jColej 
jet al.j ( |2014| ) ; [Oh et ^ ( |2015| > ; and E irby et ah] ( |2014| ). 


3.3 Dark matter cusp-core transformations 


We now turn our attention to the evolution of the under¬ 
lying dark matter halos in our simulations. In Figure^ we 
show the time evolution of the dark matter density pro¬ 
file for our M 200 = 10^ M© (left), 5 x 10® M© (middle) and 
10® M© (right) fiducial simulations. Marked on each panel is 
the 3D (ri/ 2 ) and projected {R1/2) stellar half mass radius 
of the stars (grey and black vertical dashed lines, respec¬ 
tively), and the time evolution of the dark matter density 
profile (coloured lines, as marked). The thin black dashed 
lines show a fitting function that we discuss in §3.3.3| 


As can be seen, in all cases the initial cusps are trans¬ 
formed into cores of size ~ the 3D stellar half mass radius 
ri /2 (grey vertical dashed lines). This should perhaps not 
be surprising. As demonstrated in previous works, poten¬ 
tial fluctuations due to bursty star formation (of the sort 
seen in our simulations here) produce cusp-core transfor¬ 
mations (e.g. Read fc Gilmore||2005 Pontzen & Governato 


2012 Pontzen et al.|2015|). Such fluctiations < 

scale at which stars form, which is ~ ^ 1 / 2 . 


the 


Notice that core formation proceeds much more quickly 


in the M 200 = 10® M© simulation than in the 10® M© sim¬ 
ulation. This is simply because the core is smaller and so 
the dynamical time within the core is smaller. The 10® M© 
simulation forms almost no stars after ~ 4 Gyrs - the SNe 
feedback is sufficient to shut down star formation without 
appealing to any external environmental agent. Yet core for¬ 
mation is already complete by this time. By contrast, the 
5 X 10® M© and 10® M© simulations continue to form stars 
for a Hubble time, and their cores continue to grow and hat- 
ten over this time. Interestingly, our 10^ M© simulation did 
not form enough stars for core formation to complete; the 
hnal dark matter density prohle after 14 Gyrs of evolution 
is very close to the initial NFW density prohle. However, 
we lack the resolution and physics to trust this result (only 
two star particles form). For this reason, we do not show the 
results for this simulation here. We will return to the ques¬ 
tion of whether dwarf galaxies are expected to be cusped or 
cored - or even exist at all - in 10^ M© halos, in future work. 

In the following sub-sections, we discuss the sensitivity 
of our results to our choice of numerical ‘sub-grid’ param¬ 
eters and initial conditions; and we present a convenient 
htting function that captures the size and time evolution of 
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Figure 4. Dark matter cusp-core transformations. From left to right, the panels show the time evolution of the dark matter halo density 
profiles in our 10® Mq, 5 x 10® Mq and 10^ Mq fiducial simulations, respectively (M8c28_4e5; M5e8c25_2e6; M9c22_4e6; see Table [^. 
The black lines mark the initial NFW profiles; the coloured lines show the time evolution, as marked. The vertical grey and black lines 
mark the 3D stellar half mass radius (ri/ 2 ) 2D projected stellar half mass radius {R 1 / 2 ): respectively. The thin grey dashed lines 
show the COReNFW profile fitting function (see ^3.3.3| for further details). 


this dark matter core as a function of the star formation 
time and the projected stellar half mass radius Ri/ 2 - (We 
use the projected half mass radius rather than ri /2 because 
Ri /2 is more useful observationally.) 


3.3.1 Sensitivity to our numerical ‘sub-grid’ parameters 

We have asserted so far that at the adopted minimum cell 
size of Ax = 4 pc we should become largely insensitive to our 
choice of ‘sub-grid’ numerical parameters. In this section, we 
explicitly test this by running a suite of simulations at fixed 
halo mass (M 200 = 10® M©) while varying a number of key 
parameters in relation to the fiducial simulation (see also 
Table[^. Our suite of ‘sensitivity analysis’ simulations differ 
from the fiducial run, as follows: 


• M8c28_4e6_rle3 | In this simulation, we increase the 
number of dark matter particles from N = 4 x 10^ to 4 x 10® 
while raising the star formation density threshold to p* = 
10®mu cm“®. In this way, we simultaneously test the effect 
of both resolution and sensitivity to p*. 

• M8c28_4e5_e001 | In this simulation, we decrease the 
local star formation efficiency per free-fall time (§ |2.2.1 ) from 
Cff = 10% to 1% (e.g. Krumholz Tan|2007 ); 

• M8c28_4e5_pST | In this simulation, we explicitly in¬ 
ject the momentum (pst) generated during the Sedov-Taylor 
phase for every SN explosion (see discussion in § 2.2.2), in ad¬ 
dition to the initial thermal energy. This amounts to an over¬ 
injection of momentum from SN feedback, but is nonetheless 
useful to test our sensitivity to feedback parameters. 


The resulting, almost identical, dark matter density 
profiles, after 14 Gyr of evolution and for each of the above 
parameter variations, are shown in FigureThis illustrates 
that order of magnitude changes in simulation parameters 
have no significant effect on our conclusions (see also e.g. 
Saitoh et al.||2008 Hopkins et al.||2011 ). 


In addition to the above tests, we have have rerun simu¬ 
lations where we record the gas densities during SNe events 
in order to quantify the fraction of resolved SNe events, see 


discussion in § |2.2.2[ In the case of the M5e8c25_2e6 simula¬ 
tion, we satisfy the rcooi ^ 3Ax resolution criterion for 97% 
of all SNe events. The reason for this is that pre-SN feed¬ 
back, here radiation pressure and stellar winds, as well as gas 
turbulence, clears star forming regions of dense gas before 
most SNe events occur (up to ~ 40 Myr after star forma¬ 
tion). Even though the adopted star formation threshold in 
our fiducial suite of simulations is 300 mu cm“®, only 26 out 
of almost ~ 55,000 explosions occur at m ^ 300?TiHcm“®, 
with the absolute majority (~ 99%) occurring at densities 
n < 1 mu cm“®. 


3.3.2 Sensitivity to our choice of initial conditions 


In this section, we test the robustness of our results to our 
choice of initial conditions. Recall that our simulations are 
somewhat idealised in that we start with a fully formed 
NFW dark matter halo, filled with the universal baryon 
fraction of gas. For our fiducial runs, the gas was set up as 


an NFW profile in hydrostatic equilibrium (see ^2.1). Here, 


we explore the effect of this assumption in an additional 
simulation - M5e8c25_2e6_rhocon. This is identical to our 
M 200 = 5 X 10® M© fiducial run but initialised with con¬ 
stant gas density out to r 2 oo (see Table [^. The results are 
shown in Figure The left panel compares the star for¬ 
mation rate as function of time; the right panel shows the 
resulting change in the dark matter density profile at the 
end of the simulation. 

Firstly, notice that M5e8c25_2e6_rhocon forms no stars 
for the first ~ 0.4 Gyr after which it has a substantial star- 
burst at ~ 1 Gyr. This occurs because the constant den¬ 
sity slab must first cool to form stars, unlike our fiducial 
run M5e8c25_2e6 that is initialised with a central density 
cusp. Once the constant density slab has cooled, however, 
it provides a much higher accretion rate of gas than our 
fiducial run since it is initially much cooler and denser at 
large radii. This drives a higher star formation rate for the 
whole simulation, with a correspondingly larger final stel¬ 
lar mass of M* = 2.3 x 10® M© as compared to our fidu- 
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Figure 6. Sensitivity of dark matter cusp-core transformations to the choice of initial conditions. Here we rerun our fiducial M 200 = 
5 X 10® Mq simulation (blue) using a constant initial gas density (M5e8c25_2e6_rhocon; green). The left panel compares the star formation 
rate as function of time; the right panel shows the resulting change in the dark matter density profile at the end of the simulation. 



Figure 5. Sensitivity of dark matter cusp-core transformations 
to the choice of simulation parameters. Here we re-run our fidu¬ 
cial simulation at 10® M© (red; M8c28_4e5), explicitly injecting 
the momentum (pst) generated during the Sedov-Taylor phase 
for every SN explosion (light red; M8c28_4e5_pST); decreasing the 
local star formation efficiency per free-fall time from Cff = 10% to 
1% (pink; M8c28_4e5_e001); and increasing the number of dark 
matter particles by a factor 10 while raising the star formation 
density threshold to p* = 10® mn cm“® (blue; M8c28_4e6_rle3). 
The initial NFW density profile is marked in black. The vertical 
grey and black dashed lines mark the 3D stellar half mass radius 
(ri/ 2 ) and 2D projected stellar half mass radius {Ri/ 2 )^ respec¬ 
tively. Notice that none of these changes produces a significant 
change in the central dark matter density profile; in all cases there 
is a dark matter core of size ~ ^ 1 / 2 - 


cial run that forms just M* = 0.6 x 10^ M©. This demon¬ 
strates that our results are sensitive to the cosmic gas ac¬ 
cretion history. We will study this further in future work 
where we will model isolated dwarfs in their cosmological 
context. Our main result, however, that dark matter cores 
of size ~ ri /2 form, is robust. In Figure]^ we see that both 
our fiducial run and M5e8c25_2e6_rhocon form cored dark 
matter profiles at the end of the simulation, of size ~ ri/ 2 - 
For M5e8c25_2e6_rhocon, the central dark matter density is 


lower by a factor ~ 2 as compared to our fiducial run, re¬ 
flecting the larger energy input from its larger stellar mass 
(see also ^4.1). 

Finally, we stress that M5e8c25_2e6_rhocon was delib¬ 
erately chosen to be extreme in order to test the maximal 
sensitivity to our initial conditions. In ACDM, dark matter 


halos assemble primarily through successive mergers (White 
fc Rees|1978 ). It is highly unlikely that following such an as¬ 


sembly, the gas would be arranged as a cold constant density 
slab. At leading order, the gas should trace the underlying 


dark matter as in our fiducial setup (e.g. Wetzel & Nagai 


2015). 


3.3.3 A convenient COReNFW fitting function 

In this section, we derive a convenient fitting function 
to capture the time evolution of our dark matter halos: 
the COReNFW profile. We find that the spherically aver¬ 
aged dark matter profile evolution of our dwarfs is well- 
characterised by a modified NFW functional form: 


McNFw(< t) = Mnfw(< T)f^ (16) 

where Mnfw(< r) is as in equation]^ and the function 
generates a shallower profile below a core radius Vc’. 



(17) 


where the parameter 0 < n ^ 1 controls how shallow the 
core becomes, where n = 0 corresponds to no core and n = 1 
to complete core formation. We tie the parameter n to the 
total star formation time tsF^ 

n = tanh(^) ; q = (18) 

^dyn 

where tdyn is the circular orbit time at the NFW profile scale 
radius Vs (see equation]^: 


^dyn — 27r 


GMnFw(^s 


(19) 
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Note that the ‘star formation time’ tsF here simply means 
the length of time that the simulation has run for (since 
stars form at a near continuous rate). 

The dark matter core size is set by the projected stellar 
half mass radius of the stars: 

rc = riRi/2 ( 20 ) 

which leaves us with just two tuning parameters 77 and n. We 
find Hi — 0.04 and 77 = 1.75 give a good characterisation of 
our simulation results over the full mass range (see the thin 
dashed lines in Figure]^. Notice that in the limit tsF ^ 0, 
Mcnfw ^ Mnfw and we recover the usual NFW profile. 
Similarly for r ^ rc, we return to the NFW form. This 
is advantageous as it means that our fitting function con¬ 
serves NFW profile mass; while the virial mass M200 and 
concentration parameter c (see equation take on their 
usual meanings. 

The density profile follows from the radial derivative of 
the mass profile as: 

PcNFw(^) = j Pnfw H- 4'jYr‘^r - ^nfw (21) 

Our functional form introduces two new parameters as com¬ 
pared to the NFW profile - tsF and R 1/2 - that control the 
core flattening and size, respectively. Both are observable^ 
however, allowing our derived density profile to be readily 
compared to real data. We will consider such comparisons 
in forthcoming papers. 



Figure 7. Sensitivity of dark matter cusp-core transformations 
to the halo spin parameter X'. Here we shows results for the dark 
matter density profile evolution in simulation M9c22_4e6_lam007 
(see Table the lines are as in Figure]^. This is identical to our 
fiducial M200 = 10® Mq simulation but with twice the initial dark 
matter halo spin parameter (A' = 0.07), resulting in a larger pro¬ 
jected half stellar mass radius. Notice that the coreNFW profile 
shown by the thin dashed lines qualitatively captures the evolu¬ 
tion of the dark matter density profile, but the match is poorer 
than for the fiducial simulation in Figure ^ This suggests that 
the coreNFW tuning parameters 77 and a have some weak de¬ 
pendence on halo spin A' and concentration c; we defer further 
analysis of this to future work. 


3.4 Scatter in the coreNFW tuning parameters 


Our coreNFW profile introduces two new parameters that 
are observable: the total star formation time tsF and the 
projected stellar half mass radius R 1 / 2 ] and two new pa¬ 
rameters that are tuned to fit the simulations: 77 that con¬ 
trols the size of the dark matter core via equation and 
K that controls how rapidly the core forms via equation 
One might reasonably ask whether these latter two show 
some scatter with halo concentration c, spin parameter A^ 
and/or assembly history. Our small sample of simulations is 
not sufficient to fully answer this question, but we can obtain 
some hint using our high spin simulation: M9c22_4e6_lam007 
(Table [^. This is identical to our fiducial M200 = 10® M© 
simulation but with double the cosmic mean spin parame¬ 
ter. The observable properties of this simulation are qual¬ 
itatively very similar to our fiducial simulation. However, 
perhaps unsurprisingly, it has a larger projected stellar half 
mass radius {R 1/2 = 0.81 kpc as compared to R 1/2 = 0.6 kpc 
for the fiducial simulation). This allows us to test how well 
the coreNFW form works for rather extreme changes in 
the spin parameter. (Note that A' = 0.07, as assumed in 
M9c22_4e6_lam007, is relatively rare; ~ 8% of halos will have 
A' > 0.07; e.g. Bullock et al.|200T ) The results are shown in 
Figure Notice that the coreNFW profile shown by the 
thin dashed lines qualitatively captures the evolution of the 
dark matter density profile, but the match is poorer than 
for the fiducial simulation in Figure This suggests that 
the coreNFW tuning parameters 77 and k have some weak 
dependence on halo spin A^ and concentration c; we defer 
further analysis of this to future work. 


4 DISCUSSION 

4.1 The energetics of cusp-core transformations 

Several papers in the literature have argued that cusp-core 
transformations should become energetically unfavourable 
below some critical stellar or dark matter mass (e.g. 
[Penarrubia et ah] 2012| |Garrison-Kimmel et al.||2013[ and 
see discussion in Yet here, we argue that cores always 
form if star formation proceeds for long enough. In this sec¬ 
tion, we explain why we arrive at a different result. 

Firstly, let us check that our simulations make sense 
from simple energetics arguments. The difference in gravita¬ 
tional binding energy of our COReNFW profile with respect 
to the NFW profile can be calculated as: 

AW = -^J^ ( 22 ) 

which is the amount of energy required to transform our 
initially cuspy profiles to cored profiles. We compare this 
with the available energy from SNe in Table As can be 
seen, in all cases < 1% of the available SNe energy is re¬ 
quired to unbind the cusp. Core formation does take time, 
however. If star formation is truncated then there will not 
have been enough integrated SNe energy to unbind the cusp, 
the extreme example being the case where we truncate star 
formation after a single star has formed. In this sense, it 
is absolutely correct that, for a given halo mass, there is a 
critical stellar mass at which core formation becomes ener¬ 
getically impossible. Our simulations suggest, however, that 
there is no critical halo mass at which this is the case - at 



















14 J. L Read; 0. Agertz; M. L. M. Collins 


least not down to 10^ M©. It is possible that core forma¬ 
tion ceases at ~ 10^ Mq because so few stars form, but the 
simulations we present here do not resolve this mass scale 
well enough to be able to draw strong conclusions. Below 
this mass scale, however, we expect galaxy formation to be 
extremely challenging. If ACDM is correct, we would expect 
dark matter halos below ~ 10^ M© to be almost devoid of 
stars and, therefore, to retain their pristine dark cusp. 

Let us now connect explicitly to the calculation in 


Penarrubia et al. (2012). Their equation 6 gives the super¬ 
novae energy used to transform cusps to cores as: 


AE ^ 

—— = > 8M0)eDM 

Esn {rri:,} 


(23) 


where (ttIhc) and ^ are the mean stellar mass and the frac¬ 
tion of mass in stars that go supernova, respectively; and 
cdm is the efficiency of coupling of the SNe energy to the 
dark matter (i.e. the same quantity that appears in Ta¬ 
ble [^. Both (m*) and ^ depend on the assumed IMF. 
For the [Chabrier (2003) IMF we assume here (see ^2.2.2), 
these are {m^) = 0.83 (averaged over the range 0.1 < 
m*/Mo < 100) and ^ = 0.00978. With these numbers, 
and using M*,birth (i-e. corrected for mass loss due 
to stellar evolution), we find: AE/Eq^ = [3.8,34,160] for 
M 200 = [10®, 5 X 10®, 10^] Mq, respectively. This gives excel¬ 
lent agreement with the AVF/F^sn reported in Table as it 
should. 

Penarrubia et al. ( |2012| find that core formation should 
become inefficient below M 200 ~ 10^° M©. We can now use 
equationj^to understand why we find cores at much smaller 
mass scales than this. The first potential effect is in our dif¬ 


ferent assumed IMF. If we switch to the Kroupa (2002) IMF 
assumed by Penarrubia et al. (2012), we have (m*) = 0.4 
and ^ = 0.0037. Since it is the ratio of these quantities that 
matters, this makes only a small difference, lowering the 
available SNe energy by ~ 20%. The second effect is the as¬ 
sumed core size Vc which affects the amount of SNe energy 
required for core formation (AVF). [Penarrubia et al. (2012) 
assume either rc = 1 kpc or Vc = O.lrg. For our 10^ M© halo, 
Tc = 0.83 kpc, while O.lrs = 0.1 kpc so our dark matter cores 
sit comfortably within the range assumed by |Peharrubia| 
et al.| (2012). Then there is the energy coupling efficiency 


cdm- Penarrubia et al. (2012) assume cdm = 0.4 which is 


substantially larger than we find here. However, this should 
act to make core formation much easier rather than harder. 
All of this leaves just one variable left: the stellar mass for 


a given halo mass. Penarrubia et al. (2012) assume a stellar 


mass to halo mass relation taken from the low mass ex¬ 


trapolation of the Moster et al. (2010) abundance matching 


relation. While this relation is extrapolated, it does match 
very well a recent abundance matching measurement for the 
Local Group by Brook et al. (| 2014 ) tha t reaches down to 
M* ~ 10 ® M©. The Moster et al. | 2010 ' ) relation assigns a 
stellar mass of just M* ~ 5 x 10 M© to a M 200 ~ 10 ® M© 
galaxy - substantially less than the M* ~ 3.5 x 10 ® M© that 
our simulated dwarf forms. It is clear, then, that our simu¬ 
lated dwarfs do not match the low mass extrapolation of the 


Moster et al. (2010) relation. This is the key reason why we 


are able to form cores ‘all the way down’, unlike [Peharrubiaj 
et al. (2012). Is this a problem, however? We have shown 


already that our dwarfs are remarkably realistic, giving a 


M200/M© 

10® 

5 X 10® 

10® 

N, 

86 

2721 

15,390 

M^/Mq 

2 X 10^ 

6.2 X 10^ 

3.6 X 10® 

M,.,birth/M© 

4 X 10^ 

12.6 X 10^ 

7.13 X 10® 

Nsn 

465 

14,761 

83,483 

AW/Fsn 

3.8 

33.4 

161.9 

CDM 

0.8% 

0.23% 

0.19% 


Table 2. Supernova energy coupling efficiency required to pro¬ 
duce the cusp-core transformations in our fiducial simulations. 
The rows give the dark matter halo mass of the simulated dwarf 
(M 200 ); the number of star particles formed the total stel¬ 
lar mass M*; the total stellar mass at birth M^^birth (i-c. before 
mass loss due to stellar evolution); the total number of SNe Nsn; 
the gravitational binding energy required to transform the cusp 
to a core AW, in units of the SN energy (Fsn = 10^^ ergs); and 
the coupling efficiency e = AW/{Eq-^Eq-^). Notice that for all 
simulations, the coupling efficiency is < 1%. 


good match to all known data for real isolated dwarfs ([ 3.2). 
Furthermore, a recent independent ‘tidal’ mass estimate for 
the Carina dwarf spheroidal - that has a stellar mass of 
M* = 4.8 ± 0.5 X 10® M© - puts its mass before infall onto 


the Milky Way at just M 200 = 3.6_2;3 x 10 M© (Ural et al. 


2015). Thus, Carina’s stellar mass and halo mass are in ex¬ 


cellent agreement with the simulations that we present here. 
And - just like our simulations - Carina is in strong conflict 
with the low mass extrapolation of the Moster et al. (2010) 
relation. As pointed out by Ural et al. (2015), this likely 
means that abundance matching fails below ~ 10^® M©, at 
least inside the Local Group. After all, we know that most 
of the Local Group dSphs have had their star formation 


truncated due to ram pressure and tides (e.g. McConnachie 


2012). This will act to lower M* for a given pre-infall M 200 , 


destroying the monotonic relation between and M 200 
that is an implicit assumption of the abundance matching 
machinery. We will present a much more detailed compari¬ 
son of our simulations with abundance matching constraints 
in a separate companion paper. 


4.2 Implications for the cusp-core problem: where 
to find pristine cusps 

Our results in Figure [^suggest that - if ACDM is correct - 
‘pristine’ dark matter cusps will be found either in galaxies 
that have truncated star formatiorQ or at radii r > ri/ 2 . 
This amounts to a strong prediction of the theory that can 
be distinguished from modifications to dark matter, since 
modified dark matter models will act in the same way on all 


halos irrespective of their star formation history (e.g. Rocha 
et al.|2013 ). Ohorbe et al. (2015) already noted that it takes 
time for cores to form. Similarly to our results here, they 
find that galaxies with truncated star formation should be 
more cuspy. Here we add two additional key points: (i) cores 
should always be of order the stellar half mass radius and 
thus ‘pristine’ cusps may also be found at r > ri /2 (pro¬ 
vided ri /2 < Ts); and (ii) we provide a fitting function - the 


^ Note that intermittent or ‘choked’ star formation would also 
significantly damp dark matter core formation, as might occur in 
some reionisation models (e.g. |Ricotti|200^ , or if dwarfs are not 
fully ram pressure stripped in a single pericentre passage. 
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COReNFW profile - that predicts the cusp-core evolution as 
a function of star formation time. This allows us to compare 
the model with data even for galaxies that are transitioning 
between being cusped and cored. We will consider such com¬ 
parisons between model and data in a series of forthcoming 
papers. 


4.3 Implications for the missing satellites 
problem: tidal scouring of dwarfs 

Notice that after cusp-core transformation, the central dark 
matter density of all of our dwarfs is almost the same with 
p(0) ~ 10® M© kpc“® (see Figure]^. This is important since 
such low density dwarfs will be very susceptible to tidal 


forces when falling into a larger host galaxy (e.g. Read et al. 


2006d). Such physics will be necessarily missed by simula¬ 
tions that model only the dark matter fluid since they do 
not include the physics required to model cusp-core transfor¬ 
mations ( Read et al.||200^ Peharrubia et aL]|2010 Brooks 


et al.||2013| ). 

We will consider the effect of tides in detail on these 
dwarfs in a forthcoming publication. Here, we use the tidal 


radius formula from Read et al. (2006d) to estimate at what 


mass and orbital pericentre such cored dwarfs will be com¬ 
pletely destroyed when falling into the Milky Way: 


GMgjx) GMg{x-rt) GMsjn) 

(x — rt )2 

n — 2aQQsrt 


(24) 

(25) 


where Ms{r) and Mg{x) are the satellite and host galaxy 
mass distributions, respectively; rt is the satellite tidal ra¬ 
dius; X is the distance from the host to the satellite; Q, = 
J/x^ is the angular velocity of the satellite of specific an¬ 
gular momentum J about the host; = GMs{rt)/rt is 
the angular velocity of stars within the satellite; and the 
parameter a corresponds to the orientation of orbits within 
the satellite with respect to its orbit about the host: 


{ 1 prograde 

0 radial (26) 

— 1 retrograde 

here we assume a = 1, while the specific angular momentum 
of the satellite follows from the orbital peri- and apocentre: 


j2 ^ ^ {^g(Xp) - ^g{Xa))XaXp 
Xp X(j 


(27) 


where is the gravitational potential of the host galaxy. 


Similarly to Aden et al. (2009), we crudely approximate 


the ‘Milky Way’ halo as a spherical Hernquist profile with 
mass and scale length Mh = 10^^ M©; Vh — 20kpc. For 
the satellite we assume a COReNFW profile. We define the 
satellite as being ‘destroyed’ if the tidal radius is smaller 
than the dark matter core radius: n < Vc- The core radius 
is set by the projected stellar half mass radius Vc = 1.75Ri/2 
(c.f. 13.3). We find that the relationship between R 1/2 and 
M 200 for our simulations agrees within a factor of two with 


& Kravtsov (2015b), namely: 


that found in Kravtsov (2013); Shibuya et al. (2015); Agertz 


Ri/ 2 ^ ri/2 ~ 0.015x200 


(28) 


where r 2 oo is related to the virial mass M 200 via equation 
Thus, the core radius Xc oc M^qq grows weakly with increas¬ 
ing mass. 

With the above approximations, we find that satel¬ 
lites are destroyed almost independently of mass and orbital 
apocentre up to ~ 10^° M©; all satellites with Xp ^ 30kpc 
and M 200 ^ 10^° M© will be tidally shredded rapidly after 
infall onto the Milky Way. The only way to avoid this is 
to truncate star formation and maintain a steeper central 
cusp. Indeed, the only reason why many low mass satellites 
are found around the Milky Way at all may be because their 
star formation was truncated rapidly by ram pressure strip¬ 
ping. Recent models suggest that ram pressure can remove 
almost all of the ISM of a dwarf in just one pericentric pas¬ 
sage ( Gatto et al.||2Q13 ). This may be vital for the survival 
of low mass dwarfs that orbit close to our Galaxy. Indeed, 


Peharrubia et al. (2010) noted already that if the ultra-faint 


dwarfs were cored they would not survive many orbits at 
their current locations in the Galax^j^ 

Such tidal scouring of dwarfs will dramatically reshape 
the subhalo mass function inside the Milky Way (and simi¬ 
larly for other large spiral or group environments). Surviv¬ 
ing subhalos will be those either on benign orbits that keep 
them far from the Milky Way, or those that fell in early 
and had their star formation truncated (and thus managed 
to maintain a steep central density cusp). In this context, 
it is perhaps telling that the two low mass dwarfs in the 
Milky Way with extended star formation (Garina and For¬ 
nax) have large orbital pericentres ( Lux et M^|2010 ). While 
the errors are large and proper motion measures fraught 
with difficulty, the latest data point to Fornax being on a 
cosmologically bizarre near-circular orbit; while Garina has 
a median Xp ~ 45 kpc. Typical orbits in AGDM should have 
apo-to-peri ratios closer to 1:5 for the surviving satellites 
and 1:6 for the accreted ones (e.g. Kazantzidis et ar]|2Q08 


Read et al.||2Q08| Lux et a'r]|2Q10[ ). The fact that the Milky 
Way classical dSphs have such seemingly circular orbits can 
be explained if all the dwarfs on more eccentric orbits were 
destroyed by tides. This then naturally alleviates the missing 
satellites problem, but in a manner very different to most 
other explanations in the literature. Typical solutions in¬ 
volve various prescriptions for painting stars inefficiently on 
dark matter subhalos (e.g. Koposov et al.||2009 ); in the pic¬ 


ture we present here most subhalos will simply cease to ex¬ 
ist. Such a scenario was explored recently by [Brooks et al.| 
(2013) where they show that it can indeed solve the missing 


satellites problem. 

There are three final implications of the above worth 
highlighting. Firstly, since subhalos are actually destroyed 


this will affect methods like leasing (e.g. Metcalf & Zhao 


2002 Bacon et al. 2010) or ‘satellite stream bumps’ (e.g. 
Johnston et al. 2002) that gravitationally probe the ex- 


® [Arraki et al.| ( |2014| point out that ram pressure stripping of gas 
can also lower the central stellar and dark matter density of satel¬ 
lites (see also a similar discussion in [Read &: Gilmore|20Q^ . This 
effect must be relatively small, however, if dSphs are to survive 
within the harsh tidal field of the Milky Way and/or Andromeda. 
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istence of purely dark satellites. Both should find a sub¬ 
halo mass function that is depleted (at least over the mass 
range 10 ® ^ M 200 /M 0 ^ 10 ^°) as compared to predictions 
from pure dark matter simulations. Secondly, while surviv¬ 
ing dwarfs may be on relatively benign orbits, their dark 
matter cores can still play a role in facilitating their trans- 


ticular, Kazantzidis et al. (2013) show that cored dwarfs on 


formation from dwarf irregulars to dwarf spheroidals (Mayer 
et al. 2001 1 Lokas et al. 2012 Kazantzidis et al.j2013 ). In par- 


‘benign’ orbits can still be efficiently transformed by tides. 
Thirdly, since both low and high mass satellites are effi¬ 
ciently destroyed this will also have an impact on the ‘too 
big to fail’ problem. We discuss this further, next. 

4.4 Implications for ‘too big to fail’ 

Once cusp-core transformations are taken into account, 
many massive subhalos will be tidally shredded on infall 
to the Milky Way, as discussed above. However, subhalos 
that fall in early enough will have their star formation trun¬ 
cated before core formation is complete, allowing them to 
survive. Such early infalling halos are known to be preferen¬ 
tially those that form early at high redshift (e.g. Diemand 
et al.|2005 ). Indeed, it has already been noted that the spa¬ 


tial and orbit distribution of the surviving dwarfs is more 
consistent with these early biased peaks than that of the full 


subhalo distribution (Diemand et al. 2005 Moore et al. 2006 
Lux et ar]|2010 ). This is usually attributed to dwarfs form¬ 


ing more easily before reionisation. However, several authors 
have argued that once satellites obtain some cold gas, they 
can efficiently self-shield against reionisation and continue 


to form stars (e.g. Gnedin 2000 Susa & Umemura 2004). 


It is perhaps telling that no unambiguous reionisation fea¬ 
ture in the star formation history of surviving dwarfs has 


been reported to date, despite many studies looking for one 

(e.g. |Cole et al. 2007 

Hidalgo et al. 2011 

Skillman et al. 

2014| Weisz et al.|2014 

). This suggests another explanation 


for why biased peaks are favoured as survivors. Here, we 
suggest that it is their early infall that matters. This causes 
their star formation to be truncated early by ram pressure, 
leading to the survival of their dark matter cusps. This then 
allows them to survive in the tidal field of the Milky Way. 

This still leaves us with a puzzle, however. If the clas¬ 
sical dwarfs are early infallers that had their star formation 
truncated then we seem to return once more to the ‘too- 


big-to-fail’ problem. Boylan-Kolchin et al. (2011) show that 


even early infalling subhalos are too massive to be consis¬ 
tent with the Milky Way classical dSphs. However, cusp-core 
transformations once again have an important role to play. 
Firstly, dwarfs on benign orbits like Carina and Fornax can 
form stars for a Hubble time and undergo complete core for¬ 
mation. Secondly, even dwarfs on more extreme orbits could 
undergo partial core formation. As shown in §3.3| complete 
core formation lowers the central stellar velocity dispersion 
by a factor ~ 2 inside the stellar half mass radius. In Figure 
we explore the effect of this on the ‘too big to fail’ problem. 
We plot the circular velocities within the half light radius of 
the Local Group dSphs, derived as in Collins et al. (2014). In 


the left panel, we use the measured velocity dispersions; in 
the right panel, we artificially raise the measured dispersions 
by a factor 1.5 to crudely ‘undo’ the effect of (incomplete) 
cusp-core transformations. The coloured bands in both pan¬ 


els mark dark matter halo circular speed curves taken from 
Boylan-Kolchin et al.| ( |2011| ) ; the peak circular velocity Umax 
is marked in each case. 

As can be seen in Figure compensating for cusp-core 
transformations (right panel) leads to many of the Milky 
Way (red) and Andromeda (blue) dwarfs becoming consis¬ 
tent with inhabiting Umax ~ 40km/s halos. According to 
Boylan-Kolchin et al. ( |2011| , there should be roughly 10 
dwarf galaxies around the MW(/M31) that live in such mas¬ 
sive halos. With the factor of 1.5 uniformly applied to the 
measured velocity dispersions, we end up with 13 MW dSphs 
(red) and 14 M31 dSphs (blue) that live in Umax 40 km/s 
halos. This is sufficient to completely solve the ‘too big to 
fail’ problem. Whether this works in detail will require more 
sophisticated modelling of the dwarf population in its cos¬ 
mological context that we defer to future work. Here, we 
simply point out that cored dwarfs have colder central ve¬ 
locity dispersions and that this already goes a long way to 
solving ‘too big to fail’, as has been emphasised already by 
previous authors (e.g. Read et al.|2006c Madau et al.|2014 ). 


5 CONCLUSIONS 

We have used high resolution simulations of isolated dwarf 
galaxies to study the physics of dark matter cusp-core trans¬ 
formation at the edge of galaxy formation M 200 = 10 ^ — 
10^ M 0 . We worked at a resolution (~4pc minimum cell 
size; ~25OM0 per particle) at which the impact from in¬ 
dividual supernovae explosions can be resolved, becoming 
insensitive to even large changes in our numerical ‘sub-grid’ 
parameters. Our key results are as follows: 

• Dark matter cores of size comparable to the stellar half 
mass radius ri /2 always form if star formation proceeds for 
long enough. Cores fully form in less than 4Gyrs for the 
M 200 = 10® M 0 and ~ 14 Gyrs for the 10^ M 0 dwarf. 

• We provide a convenient two parameter ‘coreNFW’ 
fitting function that captures this dark matter core growth 
as a function of star formation time and the projected stellar 


half mass radius (13.3.3). 


• We showed that our dwarf galaxies give a remarkable 
match to the stellar light profile; star formation history; 
metallicity distribution function; and star/gas kinematics of 
isolated dwarf irregular galaxies (Figures]^ and |^. In partic¬ 
ular, our results suggest that the isolated dwarf galaxy Leo 
T has a mass M 200 ~ 5 x 10® M 0 , and that we are currently 
viewing it nearly face on (at an inclination angle of ~ 20 °). 
We argue that this explains the lack of an observed rotation 
curve for Leo T. 

• We make a strong prediction that if ACDM is correct, 
then ‘pristine’ dark matter cusps will be found either in 
systems that have truncated star formation and/or at radii 
r > ri/ 2 . 

• Complete core formation lowers the projected velocity 
dispersion at ri /2 by a factor ~ 2 , which is sufficient to fully 
explain the ‘too big to fail problem’ (though we stress that 
a full solution likely also involves unmodelled environmental 
effects). 

• Cored dwarfs will be much more susceptible to tides, 
leading to a dramatic scouring of the subhalo mass function 
inside galaxies and groups. 
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Figure 8. Measured circular velocities within the half light radius of the Local Group dSphs, derived as in |Collins et al.| | |2014[ l. In the 
left panel, we use the measured velocity dispersions; in the right panel, we raise the measured dispersions by a factor 1.5 to compensate 
for (incomplete) cusp-core transformations. (Recall that cusp-core transformations can shift the central velocity dispersion of dwarfs by 
up to a factor of two; see Fi gures [^ f) and l^f).) The coloured bands in both panels mark dark matter halo circular speed curves taken 
from [Boylan-Kolchin et al.| ( [^ll| ; the peak circular velocity Umax is marked in each case. According to [Boylan-Kolchin et ar] ( |2011| >, 
there should be roughly 10 dwarf galaxies around the MW(/M31) that live in such massive halos. With the factor of 1.5 uniformly 
applied to the measured velocity dispersions, we end up with 13 MW dSphs (red) and 14 M31 dSphs (blue) that live in Umax ~ 40km/s 
halos. This is sufficient to completely solve the ‘too big to fail’ problem. 


• Our simulated dwarfs naturally lead to younger stars 
being more centrally concentrated than the older stars. This 
occurs because the older stars are collisionlessly heated sim¬ 
ilarly to the dark matter, being pushed out to larger radii. 
Such a signature is seen in the Aquarius dwarf irregular 
and is well-matched by our M 200 = 10^ M© simulation. 
Similar age-radius gradients have also been reported for a 
much larger sample of nearly dwarf irregulars by |Zhang| 


et al. (2012). They interpret the signature as an ‘outside- 


in’ shrinking of the star formation in these dwarfs. Here we 
suggest instead that it is a sign of collisionless heating caused 
by bursty star formation; the same heating that transforms 
dark matter cusps to cores. 


anonymous referee for helpful feedback that improved the 
manuscript. 


APPENDIX A: INITIAL CONDITIONS TEST 


In this appendix, we test the stability of our initial condi¬ 
tions (see ^2.1 for details of the simulation set up). In Figure 
Al we show the initial (black) and final (red) spherically av¬ 
eraged dark matter density profile for the M5e8c25_2e6 sim¬ 
ulation after 14 Gyrs of evolution, without any gas physics. 
This shows that our initial conditions are very stable. A 
small core forms due to numerical relaxation; this is of size 
~ 40 pc which is substantially smaller than any of the dark 
matter cores we discuss in this paper. 
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